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Abstract 

Two liquid state theories, the self-consistent Ornstein-Zernike equation (SCOZA) and the hierar- 
chical reference theory (HRT) are shown, by comparison with Monte Carlo simulations, to perform 
extremely well in predicting the liquid-vapour coexistence of the hard core Yukawa (HCY) fluid 
when the interaction is long range. The long range of the potential is treated in the simulations 
using both an Ewald sum and hyper spherical boundary conditions. In addition, we present an 
analytical optimised mean field theory which is exact in the limit of an infinitely long range inter- 
action. The work extends a previous one by Caccamo et al [Phys. Rev. E, 60, 5533 (1999)] for 
short range interactions. 
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I. INTRODUCTION 



Two liquid state theories are presently available which provide an accurate description 
of the liquid-vapour transition of simple systems including the critical region. A first one 
is the so-called self-consistent Ornstein-Zernike approximation (SCOZA) devised initially 
by H0ye and Stell^ 2 - and later adapted for numerical calculations^^. It is based on the 
assumption that the direct correlation function behaves, to a good approximation, at long 
range (attractive region of the potential) as the potential with a density and temperature 
dependent prefactor which is determined by imposing consistency between the energy and 
compressibility routes to the thermodynamics. The second approach put forward by Parola 
and Reatto^, incorporates renormalisation group ideas by gradually taking into account 
fluctuations of longer and longer wavelengths. It is based on the hierarchical reference 
theory (HRT) of fluids truncated at lowest order by means of an ORPA-like (optimised 
random phase approximation) ansatz (see Ref.- for review). Both predict non-classical 
critical exponents^. Thermodynamic properties and phase diagrams based on the SCOZA 
and HRT approaches have been reported for a variety of potentials (limiting ourselves to 
continuum systems) including the hard core Yukawa (HCY) potential^ 1 ^^ 2 -, the square well 
potential^ 1 ^, the Lennard- Jones (LJ) potential^ 1 ^ 1 ^ 1 ^, the Gaussian potential^, ultrasoft 
potentials^, Asakura-Oosawa pair potentials 21 ' 22 , Girifalco potentials^ 2 ^ 2 ^, binary hard 
core Yukawa mixtures 2 ^ 2 ^ 2 ^ 1 ^, etc.. Most of the comparisons with simulation data are so 
far restricted to short or medium range potentials. The aim of this article is to extend such 
comparisons to long range interactions. A convenient choice to do this is the HCY potential 
as the range can be varied from short range appropriate for modelling colloidal suspensions 
or protein solutions to medium range, where it mimics the familiar LJ potential 10 , to long 
range. For long range interactions we have considered a third approach, an optimised mean 
field theory (OMF)^, exact in the Kac limit^P, which can be worked out analytically for the 
HCY potential (see section 6). A further advantage of the HCY potential is that the repulsive 
part of the potential which in SCOZA, HRT and OMF is used as a reference potential is a 
hard sphere potential whose properties are accurately known^i^ 2 -^. Last, an approximate 
analytical solution (mean spherical approximation) is available for the HCY potential^i^ 
which drastically reduces the computational efforts of the solution of the SCOZA equations. 

For short inverse screening lengths a of the Yukawa potential, 1.8 < aa = a* < 7 (a 
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hard sphere diameter), a comparison between SCOZA, HRT and simulation results for the 
liquid- vapour coexistence curve is available inM. In this work we extend the domain of a 
from a* = 1.8 (where the HCY potential approximates the LJ potential) to a* = 0, the 
infinite long range potential. With respect to simulations, in this domain the range of the 
Yukawa potential generally exceeds the dimensions of the simulation box (for typical system 
sizes of a thousand particles) and some care is necessary to properly treat the long range of 
the potential. This has been done both by using periodic boundary conditions (b.c.) and 
performing an Ewald (EW) sum^ and by using hyperspherical b.c— 

The remainder of the paper is organised as follows: In section [TT] we define the HCY 
interaction. Section III II provides a description of the two Monte Carlo (MC) methods used 
to treat the long range of the Yukawa interaction, Ewald sum and hyperspherical method. 
After summarising the three considered theories, SCOZA, HRT and OMF in section HVl we 
compare, in section |V] the numerical results obtained with the different theories for the 
liquid-vapour coexistence curve and the thermodynamic properties along the coexistence 
boundary with the simulation results. The conclusions are presented in section IVII An 
appendix describes some properties of Yukawa charge distributions relevant for elaborating 
the OMF theory. 

II. MODEL 

In our system particles interact by the hard-core Yukawa potential 



In the limit a — > the potential gets infinitely long range and infinitely weak (Kac limit). 
We define an inverse reduced temperature (3* = l/T* = q 2 /k^Ta, where k-Q is the Boltzmann 
constant, T the temperature and q has dimension of an electric charge. The reduced density 
is p* = pa 3 . 




r < a 



r > a 



(1) 
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III. SIMULATIONS. 

Numerical simulations of fluids or plasmas involving Yukawa interactions v(r) oc 
exp(— ar)/r can be performed either in a cube of side L with periodic b.c. (C3) or 
on the surface of a four dimensional (4D) sphere of centre O and radius R (equation 
x 2 + y 2 + z 2 + u 2 = R 2 ), the hypersphere £3 for short. In both geometries the Helmholtz 
equation (A — a 2 )v(x) = — 47r5(x) can be solved analytically. 

A. Ewald sums. 

In the case of C3 the solution of the Helmholtz equation can be reexpressed in terms of 
Ewald sums with good convergence properties 3 -. The use of an Ewald potential is of course 
only necessary when the range of the Yukawa potential exceeds the linear dimensions of the 
simulation box. Molecular dynamics simulations of simple models of Yukawa plasmas were 
performed recently with this method 3 -*^. For the potential model Eq. (TjQ) the Ewald sum 
for the energy is given by 3 -^ 
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In Eq. (|2j) Tij = Tj — Ti, L is the box length, V = L 3 the volume and erfc denotes the 
complementary error function. The prime in the sum over n = (n x ,n y ,n z ), with n x ,n y ,n z 
integers, restricts it to % ^ j for n = 0. The parameter rj governs the convergence of the 
real-space and reciprocal contributions to the energy. With rj = 6.5/ L, adopted in our 
calculations, only terms with n = need to be retained in Eq. (j2J). The sum in reciprocal 
space extends over all lattice vectors k = 27rn/L with |n 2 | < 36. 



B. The hypersphere method for Yukawa interactions. 

In the 1S3 geometry the Green's function of the Helmholtz equation for the Yukawa fluid 
is obtained analytically 3 ^^: 

c , 1 sinh(u;(7r — ib)) „ 

V ^ = p • ' U { for aR > 1 > 
R smtp smh{ujip) 

1 sin(o;(7r — ib)) „ „ 

= P ■ I - ( I \ fOT aR ^ 1 ' 

where u = (|a: 2 i? 2 — l]) 1 ^ 2 . The pair potential v = —q 2 a* 2 v S3 (ip) in 1S3 is isotropic and 
depends only on the geodesic length Rip. The geodesic distance between two points Ri and 
R 2 of 1S3 is given by 

du = Ripn = -Rarccos ( ~p2~~^ ) ■ (5) 

The configurational energy of a system of N Yukawa hard spheres in £3 is therefore given 
by 

2 *2 

V(l, ...,N) = J>(^) + E t; H 3 s(^) - ^NVo . (6) 

i<j i<j 

In equation (jSj) v^ s (ipij) denotes the hard sphere potential in 5 3 , i.e. 

«| 3 S (V) =0 Rlb>a 

= 00 Rip < a . (7) 

The presence of the constant Vo in the r.h.s. of equation ([6]) accounts for the fact that 
the energy is defined up to an additive constant. The choice of Vo is a recurrent prob- 
lem of simulations in £3 (see e.g. the discussion in reference^). Here we define Vq by 
requiring that the self energies in £3 and the usual Euclidean space IR 3 coincide, i.e., 
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Vo = lim r ^ (exp(— ar)/r — v S3 (r/R)) which gives 



V = — ot + — cot(u;7r) aR < 1 
R 

= —a + — coth(c<J7r) aR > 1 
R 



(8) 



Recent Monte Carlo simulations of Yukawa plasmas have been performed with this 
methocP^. Here we apply it to a liquid of attractive Yukawa spheres. 

IV. THEORY 
A. SCOZA 

The self-consistent Ornstein-Zernike approximation in its most general formulation is a 
liquid-state theory that introduces in the relation between the direct correlation function c(r) 
and the pair potential v(r) one or more state-dependent parameters that are determined by 
imposing consistency between the different routes to thermodynamics leading to a partial 
differential equation (PDE) for the unknown parameter (s). The SCOZA originally grew 
out of the semi-analytic solution of the mean spherical approximation (MSA) for HCY 
potentials' 1 ^ and it is this SCOZA formulation that is considered in the present work. It is 
based on the Ornstein-Zernike equation supplemented with an MSA-type closure relation 

air) = for r < a 

(9) 

c(r) = c HS (r) + K(p, p)v(r) for r > o. 

Here g(r) is the pair distribution function, c HS (r) is the direct correlation function of the hard- 
core reference system given e.g. in the Waisman parameterisation 34 and K(p, 0) is the state- 
dependent parameter, that is not given in advance but is determined so that thermodynamic 
consistency is obtained between the energy and compressibility routes (/3 = l/k^T). The 
consistency requirement leads to a PDE of the form 

M^h'V' (10) 

where Xrcd = (l _ P J c(r)d 3 r) 1 is the reduced (with respect to the ideal gas) dimensionless 
isothermal compressibility given by the compressibility route and u is the excess (over ideal 
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gas) internal energy per volume provided by the energy route, i.e. u = 2iip 2 JJ° drr 2 v{r)g{r). 
Once 1/xred is expressed in terms of the excess internal energy the PDE eq. (1101) can be 
transformed into a PDE for u 

For the general case of an arbitrary attractive tail of the pair interaction the determination 
of B(p,u) must be done fully numerically^ 4 - and requires an enormous amount of compu- 
tational cost. However, this procedure is not necessary for the fluid considered here: for 
certain pair potentials, like the hard-core Yukawa potential or hard-core potentials with 
linear combinations of Yukawa and exponential tails - so called Sogami-Ise fluids, the coef- 
ficient B(p, u) can be obtained semi- analytically due to the availability of the semi-analytic 
MSA solution for these potentials. Details of the determination of B(p, u) can be found in 
Ref.-,— and 4 ^. While the former formulation of SCOZA in Ref.- is based on the Laplace 
transformation route for solving the underlying MSA for hard-core Yukawa fluids, we have 
used in the present work the reformulation described in Refs.— and 4 ^ which is based on the 
Wertheim-Baxter factorisation method for the solution of the MSA. While both solution 
methods are equivalent for the potential considered here, the Wertheim-Baxter factorisation 
method is more flexible and allows a broader applicability of the SCOZA 

The PDE is a quasilinear diffusion equation that has been solved numerically via an 
implicit finite-difference algorithm in the region (J3*,p*) G [0, /3/] x [0,1] with suitable ini- 
tial and boundary conditions (see Ref.- for more details). We have chosen a density and 
temperature grid spacing of Ap* = 10~ 4 and A/3* = 2 ■ 10~ 4 . Up to now SCOZA has been 
applied to a number of discrete^ 4 - and continuum systems&2£ and the results showed - when 
compared with computer simulations - that the theory yields very accurate predictions for 
the thermodynamic and structural properties throughout the temperature-density plane, 
even near the phase coexistence and in the critical region. Indeed, SCOZA is one of very 
few liquid state theories that do not show serious problems in the critical region and even 
exhibit some form of scaling with non-classical, partly Ising-like critical exponents 8 . 

B. HRT 

In HRT the pair potential v(r) is divided into a repulsive reference potential part vr(t) 
and a longer-ranged attractive part w(r). In the present case the reference system is the 
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hard-sphere fluid, whose thermodynamics and correlations are accurately described by the 
Carnahan-Starling equation of stated and the Verlet-Weis parameterisation of the two-body 
radial distribution function™ The HRT differs from the conventional liquid state approaches 
in the way the attractive perturbation is dealt with: in order to take into account the effect 
of fluctuations, w(r) is introduced gradually via a sequence of Q systems whose interaction 
potential vq{t) results from the sum of vr(t) and a perturbation wq(t). The latter term 
is defined by its Fourier transform Wq(k) which coincides with w(k) for k > Q, and is 
identically vanishing for k < Q. As Q evolves from Q = oo to Q = 0, the interaction 
v Q( r ) = v n{ r ) + w q{ t ) goes from the reference part Vr(t) to the full potential v(r). In other 
words the parameter Q plays the role of an infrared cutoff whose effect consists in depressing 
fluctuations on length scales larger than 1/Q: the critical fluctuations are recovered only 
in the limit Q —>■ 0. The corresponding evolution of thermodynamics and correlations of 
increasing order can be determined via perturbation theory and is described by an exact 
hierarchy of integro-differential equations. As mentioned above close to a critical point and 
at large length scales, this hierarchy becomes equivalent to a formulation of the momentum- 
space renormalisation group^. Unlike current liquid-state theories HRT preserves the correct 
convexity of the free energy also below the critical temperature, producing flat isotherms in 
the coexistence region which correspond to infinite compressibility and constant chemical 
potential. The first equation of the hierarchy gives the evolution of the Helmholtz free energy 
Aq of the partially interacting system in terms of its two-body direct correlation function 
in momentum space cg(fc) and the full perturbation w(k): 



where = —w{k)/k-QT. The quantities Aq, Cq(/e) are linked to Aq and c<g(fc) by the 
relations 



where V is the volume and p the number density. These modified free energy and direct 
correlation function have been introduced in order to remove the discontinuities which appear 
in Aq and CQ(k) at Q = and k = Q respectively as a consequence of the discontinuity of 
WQ{k) at k = Q. Indeed, they represent the free energy and direct correlation function of 




(12) 
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the fully interacting fluid as given by a treatment such that the Fourier components of the 
interaction with wavelengths larger than 1/Q are not really disregarded, but approximately 
taken into account by mean-field theory. In the Q — > limit the modified quantities coincide 
with the physical ones, once the fluctuations have been fully included. For Q = oo the 
definitions in Eqs. (1131) and (|14p reduce, instead, to the well-known mean-field expressions 
of the free energy and direct correlation function. 

A simple approximation scheme consists in truncating the hierarchy at the first equation, 
supplementing it with a suitable closure relation based on some ansatz for Cq(&). As in 
previous applications, our form of Cq{k) has been inspired by liquid-state theories, 

c Q (k) = c HS (k) + a q + g Q (k) , (15) 

where chs(^) is the Fourier transform of the direct correlation function of the hard-sphere 
fluid, and Aq, Gq{k) are a priori unknown functions of the thermodynamic state and of Q. 
The function Qq(k) is determined by the core condition, i.e., the requirement that the radial 
distribution function gg(r) is vanishing for every Q whenever the interparticle separation is 
less than the hard-sphere diameter a. Via Eq. ([15!) we can write the core condition in terms 
of an integral equation for cg(fc). Aq is adjusted so that Cq{U) satisfies the compressibility 
rule. For each Q-system this constraint gives the reduced compressibility of the fluid as the 
structure factor evaluated at zero wavevector, and can be expressed in terms of the modified 
quantities Aq, Cq(&) as 

Cq(^ = 0) = ^. (16) 

This equation can be regarded as a consistency condition between the compressibility route 
and a route in which the thermodynamics is determined from the Helmholtz free energy as 
obtained from Eq. (fl2l) . Here the compressibility rule (fl6l) plays a fundamental role. In fact, 
when Aq in Eq.( fT5|) is determined via Eq.( fl6i) and the resulting expression for Cq(/c) is 
used in Eq. (1121) . one obtains a partial differential equation for Aq which reads 



3Aq = _H_, 

dQ 4tt 2 



HQ) 



(17) 



A^(Q)+^(Q) 
where we have set Aq = d 2 A Q /dp 2 } ip{k) = $(/c)/$(0) and 

= c HS (£0 + e Q (k) - [chs(O) + G Q (o)Mk) . (18) 
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Eg. ([171) is integrated numerically from Q = oo down to Q = 0. At each integration step, 
Gci{k) is determined by the core condition gq(r) — 0, < r < a. This condition acts as an 
auxiliary equation for determining ijj(k) via Eq. fllSI) . 

For this specific work we considered the density interval p* G [0, 1] to numerically solve 
the differential equation for the free energy. The spacing of the density grid is Ap* = 10~ 3 , 
consequently the error in the determination of the coexistence densities is of the same order 
of magnitude. Concerning the boundary conditions: for low density the system behaves as 
an ideal gas while at p* = 1 we assumed the validity of the standard ORPA approximation. 

Two remarks are worth mentioning: First, Eq. ( |T5l) assumes that the fluid direct cor- 
relation function depends linearly on the perturbation $(A;). This ansatz is especially ap- 
propriate when the perturbation range is much longer than that of the reference part of the 
interaction, i.e. longer than the particle size. The second remark relates to the implemen- 
tation of the core condition: the inverse Fourier transform of the function Qq{k) has been 
expanded in a series of Legendre polynomials for r < a and the series has been truncated 
after a finite number of terms (typically five). The evolution equations for the expansion 
coefficients were then determined by differentiating with respect to Q the integral equation 
for cq(/c), i.e. the core condition, and subsequently projecting it on the polynomials used in 
the expansion. However, the resulting equations are coupled to the evolution equation [12] 
for the free energy, and this makes them difficult to handle. Therefore, as described in detail 
in [21, 22], in the derivative of cq(/c) with respect to Q the long-wavelength contributions 
containing the isothermal compressibility of the Q-system were disregarded. Physically this 
amounts to a decoupling of the short- and the long-range evolutions of the correlations. 
This approximation, as well as the previous one, are fully justified for long-range pertur- 
bations, where the short- and the long-range parts of the correlations are mainly affected 
by the reference and the perturbation terms respectively, but become more problematic for 
short-range interactions where the interplay between excluded-volume and cohesion effects 
is much stronger. However, this effect does not concern the present study. 

C. Optimised mean field theory 

We consider, quite generally, a fluid of identical hard spheres of diameter a with additional 
isotropic pair interaction u(ry). Since v(r) is an arbitrary function of r for r < a, one can 
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assume that v(r) has been regularised in the core in such a way that its Fourier transform 
v(k) is a well behaved function of k and that v(0) is a finite quantity. 

At equilibrium at inverse temperature f3 and chemical potential p (grand canonical (GC) 
ensemble) the pressure p of the fluid is a convex function of v = (3p (at fixed f3)^2- even 
before the thermodynamic limit has been taken^^. As a consequence (3p{v) is continuous, 
its derivatives with respect to v exist and it admits a Legendre transform (3f (p) with respect 
to v defined as 

Pf(p) = sup {pu-(3p{u)} . (19) 

v 

The GC free energy (3f (p) is then a convex function of the density and its Legendre transform 
with respect to p coincides with fipiy). 

It was shown in ref.— that in the case of an attractive interaction (i.e. w(k) = —/3v(k) > 
for all q) we have the inequality 

/3/GM < PfMP>p) 

PfuFiP, P) = Pfns(p) - \p 2 w(0) + \pu>{0) • (20) 

Here the subscript mean field (MF) emphasises that /mf is the van der Waals free energy 4 ^. 
Note that /3/mf [p] is a priori non convex, notably in the two-phase region. 

^From here on we specialise to the Yukawa interaction w(r) = f3q 2 a* 2 y(r) with y(r) = 
exp(— ar)/r. Since y(0) = oo, equation (1201) does not appear to be very useful. However, it 
follows from appendix A that we can replace y(r) by the regular function W r {r}IQ 2 r where 
W T (r) is the interaction energy of two spherical distributions r(r) of Yukawa charges of 
effective charge Q T . The interaction is regularised in the core but remains unchanged for 
r > a. Since W T (k) > 0, Eq. ( |20l) is still valid and will give a non trivial upper bound for 
the free energy. Similar to what has been done in reference^ for the Coulomb interaction 
we seek a best upper bound for Pf((3, p) by minimising the quantity 

W T (0)-pW T (0) 
X ~ Ql • (21) 

To this end we consider variations of X with respect to the charge distribution r(r) defined 
by Eq. ( 1A1I) and work out the stationary condition (a = a/2) 

— — - = for r < a . (22) 
ot{t) 
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First, it follows from equation (1A5|) that 

SQ T sinh(ar) 



(23) 



6t(t) ar 

Second, as a consequence of the convolution relations V T = r * y and W T = r*V T = r*y*r 
one has 

and thus 

Finally, since W T (0) = 47rcv~ 2 r 2 (0), we have 

8W T (0) _ 87rf(0) 
8t(t) a 2 

Collecting results (1231. (125]) and (EE} one gets 

<5X 2 smh(ar 



cV(r) Q 3 ar 



W T (0)-pW T (0) 



(26) 



H-if^Cr)-?^! forr<^. (27) 



Q 2 r\ ) 

Note that we can impose r(0) = 1 since, if r(r) is solution of the stationary condition (|22l) . 
then Ar(r), where A 7^ is an arbitrary constant, is also a solution (Q T is then multiplied 
by A but W T /Q 2 is left unchanged). Therefore Eq. (1221) can be recast in the form 



VAr) = + ± + !^ (»V(0) - ^(0)) . (28) 
cr Q r ar V / 



Of course, Eq. (j2~8"]l is valid only for r < a. The first contribution in the r.h.s of Eq. ( 1281 
is clearly identified with the potential created by a uniform density of Yukawa charges of 
density p while the second contribution stems from a charge discontinuity at r = a. The 
solution is therefore of the form r(r) = pO(cf) +X5(r — a), where K is a constant determined 
by the condition r(0) = 1, yielding 

r(r) = p6(ef) + ^— ^5(r - a) (29) 

(77 = 7rpcr 3 /6 packing fraction). Furthermore, from Eq. (]A5j) it follows that 

a?xcosh(acx) — sinh(a(x) sinh(acx) 
Qr = 3ri — 1-^ = — 30 

(aa) A aa 
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so that the potential V T is seen to be the superposition of two potentials created by spherical 
surface and volume distributions, the expressions of which are given in appendix A by Eqs. 
(IA13p and fjA15j) . respectively. One finds 

V T (r) = —!- — (rja* 2 + 12rj + 6r]a* - a* 2 ) * — '- . 31 

a* 2 a* 2 v 1 ar 

One can check that V T indeed verifies Eq. ( I28l) . Finally, the optimised upper bound for the 
free energy can then be shown to be equal to 

< /3/omf(/3,p)=/3/hs(p)+A mf(/3,p) 
6f3ria* s e~ a * 

AoMF — X 

7T 

r/a* 2 + 6r]a* + 12r] - a* 2 
(12/7 — 6r/a* — a* 2 + a* 2 //) — e~ a *(a* 2 — a* 2 ?? — 6r/a* — 12r/) 
-TV M-5-5,+fl a , 2 + 0(Q , 3) 
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An optimised mean field (OMF) theory for the fluid can now be devised by considering 
the Landau function^ 

£(P, Pi v) = up- PfoMF(P, p) ■ (34) 

At given (3 and v, the density p is the one which minimises the Landau function £(/3,p, v). 
In the limit a* — > 0, Aqmf = —2np 2 f3 and one thus revovers the free energy of the Kac 
model for a = 0. 



V. RESULTS 



A comparison between the three theoretical approaches and simulation data for the liquid- 
vapour coexistence curve is made in Figs. IVII - IVII for a* = 1.8, 1.0,0.5 and 0.1. Values for 
the densities, internal energies, chemical potentials and Helmholtz free energies along the 
coexistence curve are presented in Tables IVII - IVII These data may be valuable for reference 
to future approaches. 

The simulation results have been obtained with the Gibbs ensemble Monte Carlo (GEMC) 
method^ 1 ^ using mostly a total of 1000 particles (EW b.c.) and 2000 particles (S3 b.c). 
The number of cycles generated after equilibration varied from 7-40 xlO 6 according to 
temperature, each cycle comprising with equal probabilities, displacement of the iV particles, 
volume change, or exchange of 500 (EW b.c.) or 1000 particles (63 b.c). Error bars (two 
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standard deviations) were calculated by averaging densities, energies or histograms over 
blocks of 1 million cycles. Expectedly, uncertainties are largest near the critical temperature 
where fluctuations of density of the coexisting phases are important and hence convergence 
slow. Finite size effects additionally affect the critical behaviour, but also influenced to 
some extent the determination of the coexistence density on the liquid side at least for the 
smallest value a* = 0.1. Thus, for the S3 b.c, an increase of the number of particles from 
1000 to 2000 typically lowered p\ by about 3%. The opposite trend was observed with the 
EW b.c. where p\ increased by 1 — 2% by an increase of N from 1000 to 1872 (cf. table [VTl) . 
It is likely that, for a* = 0.1, small finite size effects are still present in both sets of results. 

The agreement between HRT, SCOZA and MC simulations turns out to be extremely 
good over the range of a values considered. Small differences occur only at the highest 
and lowest values of a. At a* = 1.8 the critical temperature of HRT is slightly below that 
predicted by SCOZA. Previous work 10 showed, however, that by further increasing the value 
of a the HRT critical temperature eventually gets higher than that of SCOZA. At a* — 0.1 
the coexistence curves, including the critical regions, of SCOZA and OMF are very close 
(cf. Fig. IVip ; as the OMF theory is expected to be nearly exact at this small value of a 
it is very likely that SCOZA becomes exact also in the Kac limit, though no formal proof 
seems to be available. This is further confirmed by the close agreement between critical 
parameters of SCOZA and OMF at a* = 10 -5 (see table IVl] ). The coexistence curve of HRT 
is found to be slightly narrower than that of SCOZA though in the near critical region the 
agreement between the two theories is quite good. As pointed out in ref.— the HRT flow 
cannot be solved conveniently in the Kac limit. Inspection of Fig. IVII shows that the OMF 
theory rapidly deteriorates with increase of a. Figure IVII compares HRT, SCOZA, OMF 
and MC simulation results for the internal energy along the coexistence curve at a* = 0.5. 
The agreement is similar to that for the densities. 

The variation with a of the critical temperature and density for SCOZA and HRT, 
including data of ref.— for a* > 1.8, is given in table IVII No attempt has been made to 
determine critical data in the simulations by fitting data away from the critical point by some 
power law expression. Because of the finite size effects in the critical region uncertainties of 
the results would be large and of no use to validate the theoretical predictions. In the GEMC 
method the order parameter M = pi~p g is believed to have MF behaviour (critical exponent 
(3 = 1/2) in the critical region once the correlation length exceeds the linear dimension L 
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of the simulation box . In fact we find that for a* < 1.0, M 2 varies linearly over the whole 
temperature region considered yielding an effective exponent /3 e // = 1/2). Such an analysis 
was not conclusive for a* = 1.8. 

VI. CONCLUSION 

We have compared the predictions for the liquid- vapour coexistence curve of a long-range 
Yukawa fluid obtained from advanced theoretical approaches with Gibbs ensemble Monte 
Carlo simulations. Concerning the simulations care has to be taken when the range of the 
potential exceeds the length of the simulation box. This was done by performing an Ewald 
sum in the case of a cubic simulation box and by using hyperspherical boundary conditions. 
The theoretical approaches that we applied comprised the self-consistent Ornstein-Zernike 
approximation (SCOZA), the hierarchical reference theory (HRT), as well as an optimised 
mean field theory (OMF). While the OMF yields the exact result in the limit of infinite 
range of the potential, it deteriorates with decreasing interaction range. On the other hand, 
HRT and SCOZA turn out to be in perfect agreement with simulation results over the whole 
interaction range considered. This study complements a recent comparison— for the case of 
intermediate and short range Yukawa fluid. 
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Table 1. Coexistence densities, internal energy U, excess chemical potential p, and Helmholtz free energy A at coexistence for a* — 1.8. A re ? is 

the Carnahan-Starling HS free energyH. 



T* p* pi ((7/iVfeT)„ {U/Nk.T) l 11/kT (Aa)J (Aa)f 

SCOZA 0.174 0.464 -1.055 -2.650 -2.737 

0.635 HRT 0.181 0.458 -1.106 -2.607 0.963 2.521 

MC-BW 0.175(10) 0.430(15) -1.06(4) -2.43(3) -2.74(1) 

SCOZA 0.161 0.485 -0.982 -2.770 -2.765 

0.63 HRT 0.166 0.477 -1.030 -2.735 0.890 2.657 

MC-BW 0.171(8) 0.470(15) -1.05(2) -2.67(2) -2.76(1) 

SCOZA 0.138 0.515 -0.862 -2.996 -2.825 

0.62 HRT 0.142 0.509 -0.901 -2.967 0.774 2.900 

MC-BW 0.149(5) 0.501(4) -0.94(2) -2.92(2) -2.83(1) 

SCOZA 0.105 0.565 -0.686 -3.407 -2.953 

0.60 HRT 0.108 0.560 -0.724 -3.386 0.609 3.333 

MC-EW 0.109(2) 0.553(4) -0.73(1) -3.34(2) -2.96(1) 

SCOZA 0.082 0.606 -0.556 -3.803 -3.094 

0.58 HRT 0.083 0.602 -0.578 -3.787 0.486 3.741 

MC-BW 0.083(2) 0.594(4) -0.58(1) -3.74(2) -3.10(1) 

SCOZA 0.063 0.643 -0.452 -4.204 -3.251 

0.56 HRT 0.065 0.640 -0.490 -4.192 0.396 4.155 

MC-BW 0.063(2) 0.631(4) -0.46(1) -4.13(3) -3.25(3) 

SCOZA 0.049 0.676 -0.367 -4.614 -3.424 

0.54 HRT 0.050 0.674 -0.390 -4.606 0.317 4.573 

MC-BW 0.050(2) 0.670(4) -0.39(1) -4.57(2) -3.47(8) 

SCOZA 0.0284 0.738 -0.235 -5.503 -3.834 

0.50 HRT 0.029 0.736 -0.258 -5.500 0.201 5.469 

MC-EW 0.029(2) 0.733(4) -0.25(1) -5.47(3) -3.76(10) 
a (Aa)„ = 7J j pr{ A-A-'f) v 

b ( Aa )l = NTPr( A - Ar °')l 

APPENDIX A: SOME PROPERTIES OF YUKAWA CHARGE DISTRIBU- 
TIONS 

The electrostatics of Yukawa charge distributions is similar to, but not identical with the 
electrostatics of usual Coulomb charge distributions. In this appendix we extend previous 
results obtained for special forms of spherical distributions of Yukawa charges^ 1 ^ 1 ^ 1 ^ 1 ^ to 
a general spherical distribution r(r) which satisfies: 

r(r) = for r>lf=a/2, 

r(r) > for < r < a = a/2 (Al) 

We denote by y(r) the Yukawa potential created by a unit point charge. We have y(r) = 
exp(— ar)/r and for its Fourier transform y(k) = 4n/(a 2 + k 2 ). The Yukawa potential 
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Tabic 2. Same as tabic 1 but for a* — 1.0 



T* pi pf (U/NkT) v (U/NkT) l fi/kT (Aa)^ (Aa)f 





SCOZA 


0. 


174 


0. 


396 


-0 


.987 


-2 


209 


-2.773 






0.905 


HRT 


0, 


176 





.393 


-1 


.000 


-2 


.195 




0.955 


2.169 




MC-BW 





.183(2) 





395(10) 


-1 


.025(50) 


-2 


19(5) 


-2.77(1) 








SCOZA 


0, 


.135 


o. 


450 


-0 


.788 


-2 


.583 


-2.864 






0.88 


HRT 





.136 


o 


448 


-0 


.767 


-2 


477 




0.758 


2.557 




MC-BW 





.135(2) 





439(5) 


-0 


.78(2) 


-2 


52(2) 


-2.87(1) 








SCOZA 


{) 


112 


(1 


484 


-0 


.676 


_2 


852 


-2.944 






0.86 


HRT 


o. 


114 


o. 


482 


_o 


.687 


_2 


.842 




0.650 


2.825 




MC-BW 





.112(2) 





.474(4) 


-0 


.67(1) 


-2 


79(2) 


-2.95(1) 








SCOZA 


{) 


.095 


(J 


514 


-0 


.586 


-3 


107 


-3.028 






0.84 


HRT 


{) 


.096 


(J 


512 


_o 


593 


_3 


098 




0.561 


3.083 




MC-BW 





.096(2) 





508(5) 


-0 


.59(2) 


-3 


.07(3) 


-3.04(2) 








SCOZA 


o. 


0806 


o. 


541 


-0 


.511 


-3. 


360 


-3.119 






0.82 


HRT 


o 


()81 


o 


540 


-0 


513 


-3 


.353 




0.485 


3.341 




MC-BW 





.085(3) 





543(6) 


-0 


.54(2) 


-3 


.37(3) 


-3.11(2) 








SCOZA 


o 


0685 


o 


567 


-0 


.440 


-3 


.61 5 


-3.216 






0.80 


HRT 


o. 


.069 


o 


566 


-0 


.451 


-3. 


.610 




0.424 


3.600 




MC-BW 





.070(2) 





564(4) 


-0 


.447(15) 


-3 


60(3) 


-3.21(1) 








SCOZA 


o. 


.058 


o 


.591 


-0 


.389 


-3 


.873 


-3.321 






0.78 


HRT 


o 


()59 


(J 


.590 


-0 


.401 


-3 


.868 




0.372 


3.859 




MC-BW 





.062(2) 





594(6) 


-0 


.420(15) 


-3. 


89(3) 


-3.30(2) 








SCOZA 





049 





614 


-0 


339 


-4 


140 


-3.434 






0.76 


HRT 





050 





.613 


-0 


.353 


-4. 


.133 




0.324 


4.124 




MC-BW 





.050(2) 





.612(3) 


-0 


.344(5) 


-4 


.120(15) 


-3.43(1) 








SCOZA 





041 





636 


-0 


.294 


-4. 


413 


-3.557 






0.74 


HRT 





042 





.635 


-0 


.304 


-4 


407 




0.280 


4.398 




MC-BW 





.042(2) 





634(4) 


-0 


.30(1) 


-4. 


440(25) 


-3.55(1) 
















a (Aa)„ 




1 (A 
NkT y 


- A 











V T = t * y created by this distribution at point R is given by the convolution of r and y 



V T (R) = [ dr 47rr 2 T(r)I(R,r) 
Jo 



I(R,r) = /^ eXP( ^' R r y r)l) ; (A2) 
where dfl r denotes the infinitesimal spherical angle about vector r. The integral I(R,r) 



rcada 37 ' 54 ' 55,56,57 



T , smh(ar) exp(-<xR) 

I(R,r) = — - — for R>r, 

ar R 

sinh(a.R) exp(-ar) „ . 

= ^— —— — for r > R . (A3) 

aR r 
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Tabic 3. Same as tabic 1 but for a.* — 0.5 



T* 




Pv 


pi 


(U/NkT) v 


(U/NkT) l 


fi/kT 




(Aa)J> 




SCOZA 


0.164 


0.367 


-0.923 


-2.063 


-2.824 






1.045 


HRT 


0.163 


0.371 


-0.916 


-2.086 




0.908 


2.082 




MC-EW 


0.169(8) 


0.355(8) 


-0.93(3) 


-1.99(2) 


-2.82(1) 








MC-S3 


0.182(5) 


0.384(3) 


-1.02(3) 


-2.15(3) 










SCOZA 


0.157 


0.377 


-0.885 


-2.131 


-2.839 






1.04 


HRT 


0.156 


0.381 


-0.881 


-2.153 




0.873 


2.149 




MC-EW 


0.160(8) 


0.370(8) 


-0.89(4) 


-2.05(4) 


-2.84(1) 








MC-S3 


0.165(6) 


0.384(6) 


-0.93(3) 


-2.172(22) 










SCOZA 


0.133 


0.412 


-0.764 


-2.377 


-2.899 






1.02 


HRT 


0.132 


0.415 


-0.760 


-2.394 




0.752 


2.390 




MC-EW 


0.130(6) 


0.400(6) 


-0.75(2) 


-2.31(4) 


-2.90(1) 








MC-S3 


0.133(2) 


0.415(3) 


-0.77(1) 


-2.39(2) 










SCOZA 


0.114 


0.442 


-0.670 


-2.602 


-2.962 






1.00 


HRT 


0.114 


0.445 


-0.670 


-2.620 




0.662 


2.617 




MC-EW 


0.118(4) 


0.439(4) 


-0.689(15) 


-2.58(2) 


-2.96(1) 








MC-S3 


0.111(1) 


0.441(2) 


-0.655(50) 


-2.60(1) 










SCOZA 


0.099 


0.469 


-0.594 


-2.817 


-3.030 






0.98 


HRT 


0.099 


0.470 


-0.594 


-2.826 




0.587 


2.824 




MC-EW 


0.103(4) 


0.467(4) 


-0.62(1) 


-2.74(2) 


-3.02(1) 








MC-S3 


0.095(1) 


0.466(2) 


-0.573(54) 


-2.800(15) 










SCOZA 


0.086 


0.493 


-0.528 


-3.028 


-3.101 






0.96 


HRT 


0.086 


0.494 


-0.527 


-3.035 




0.520 


3.032 




MC-EW 


0.085(4) 


0.485(4) 


-0.52(1) 


-2.99(1) 


-3.11(1) 








MC-S3 


0.0846(11) 


0.494(3) 


-0.521(7) 


-3.037(15) 










SCOZA 


0.075 


0.516 


-0.471 


-3.239 


-3.178 






0.94 


HRT 


0.076 


0.516 


-0.476 


-3.240 




0.469 


3.238 




MC-EW 


0.076(2) 


0.512(2) 


-0.468(5) 


-3.21(1) 


-3.17(2) 








MC-S3 


0.072(1) 


0.509(2) 


-0.442(4) 


-3.20(1) 










SCOZA 


0.0656 


0.538 


-0.420 


-3.451 


-3.259 






0.92 


HRT 


0.066 


0.537 


-0.422 


-3.447 




0.417 


3.446 




MC-EW 


0.068(2) 


0.539(3) 


-0.43(1) 


-3.46(2) 


-3.27(3) 








MC-S3 


0.0678(13) 


0.547(3) 


-0.437(8) 


-3.52(2) 










SCOZA 


0.0572 


0.559 


-0.374 


-3.667 


-3.346 






0.90 


HRT 


0.058 


0.557 


-0.380 


-3.658 




0.374 


3.656 




MC-EW 


0.058(2) 


0.556(3) 


-0.375(10) 


-3.64(2) 


-3.34(1) 








MC-S3 


0.058(1) 


0.566(3) 


-0.382(7) 


-3.72(2) 














S (A«)„ 
"(H 


~ NkT — 


^ re/ )„ 









The potential V T (R) for R < a is not particularly useful but outside the sphere (i.e. for 
R>a) one infers from equations (1A2I) and (1A3jl the remarkable result 



V r (R) = Q r e -^^- , (A4) 
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Tabic 4. Same as tabic 1 but for a* — 0.1 



T" 


pZ 


Pi 


(U/NkT) v 


(U/NkT), 




(Ao)5 (Aa)J> 


SCOZA 


0.158 


0.357 


-0.900 


-2.030 


-2.863 




1.10 HRT 


0.161 


0.353 


-0.916 


-2.010 




-0.916 -2.010 


MC-BW 


0.17(1) 


0.340(15) 


-0.96(2) 


-1.96(5) 


-2.86(2) 




MC-S3 


0.167(5) 


0.361(5) 


-0.93(3) 


-2.05(2) 







SCOZA 


0. 


145 


0.375 


-0.834 


-2 


.155 


-2.891 


HRT 


0. 


148 


0.371 


-0.850 


-2 


.132 




MC-BW 


0. 


145(10) 


0.360(15) 


-0.82(1) 


-2 


.05(1) 


-2.89(1) 


MC-BW1 


0. 


145(4) 


0.362(3) 


-0.82(2) 


-2 


.07(2) 


-2.895(5) 


MC-S3 


0. 


160(6) 


0.392(3) 


-0.92(2) 


-2 


.25(2) 





SCOZA 


0. 


134 


0.391 


-0.778 


-2 


.270 


1.08 HRT 


0. 


138 


0.387 


-0.800 


-2 


.244 


MC-BW 


0. 


136(5) 


0.380(5) 


-0.79(2) 


-2 


.21(2) 


MC-S3 


0. 


1384(24) 


0.395(3) 


-0.8018(15) 


-2 


.29(2) 



-0.800 -2.244 

-2.92(2) 



SCOZA 


0. 


124 


0.406 


-0 


730 


-2 


.379 


-2.947 


HRT 


0. 


128 


0.401 


-0 


749 


-2 


.347 




MC-EW 


0. 


126(2) 


0.393(3) 


-0 


73(1) 


-2 


.32(2) 


-2.95(1) 


MC-BW1 


0. 


127(2) 


0.399(3) 


-0 


74(1) 


-2 


.34(2) 


-2.95(1) 


MC-S3 


0. 


1260(15) 


0.4103(21) 


-0 


737(8) 


-2 


.402(13) 





-0.749 -2.347 



SCOZA 


0.116 


0.420 


-0.686 


-2.484 


-2.977 




1.06 HRT 


0.120 


0.415 


-0.708 


-2.452 




-0.708 -2.452 


MC-BW 


0.121(4) 


0.416(4) 


-0.716(15) 


-2.46(2) 


-2.97(1) 




MC-S3 


0.118(2) 


0.423(3) 


-0.697(11) 


-2.544(15) 







SCOZA 


0. 


1082 


0.434 


-0.647 


-2 


.587 


-3.007 


HRT 


0. 


112 


0.427 


-0.667 


-2 


.547 




MC-BW 


0. 


109(2) 


0.425(3) 


-0.65(1) 


-2 


.54(2) 


-3.000(15) 


MC-BW1 


0. 


113(2) 


0.433(3) 


-0.68(1) 


-2 


.59(2) 


-3.00(1) 


MC-S3 


0. 


111(2) 


0.442(2) 


-0.660(7) 


-2 


.63(2) 





SCOZA 


0. 


101 





446 


-0.610 


-2. 


689 


1.04 HRT 


0. 


107 





439 


-0.638 


-2. 


644 


MC-EW 


0. 


102(1) 





439(2) 


-0.641(5) 


-2. 


64(1) 


MC-S3 


0. 


1047(13) 





456(3) 


-0.630(9) 


-2. 


75(2) 



-0.638 -2.644 

-3.04(1) 



SCOZA 


0.089 


0.470 


-0.546 


-2.887 


-3.103 




1.02 HRT 


0.093 


0.462 


-0.571 


-2.837 




-0.570 -2.837 


MC-EW 


0.091(1) 


0.466(2) 


-0.555(5) 


-2.86(1) 


-3.10(1) 




MC-S3 


0.091(1) 


0.479(2) 


-0.558(6) 


-2.94(1) 







SCOZA 


0.078 


0.493 


-0.489 


-3.085 


-3 


.173 


HRT 


0.083 


0.483 


-0.519 


-3.026 




-0.519 -3.026 


1.00 MC-BW 


0.077(3) 


0.483(5) 


-0.48(1) 


-3.028(25) 


-3 


.17(1) 


MC-EW1 


0.078(2) 


0.487(5) 


-0.484(14) 


-3.05(3) 


-3 


.185(15) 


MC-S3 


0.079(1) 


0.499(2) 


-0.492(6) 


-3.126(13) 






SCOZA 


0.0572 


0.559 


-0.374 


-3.667 


-3 


.346 


0.98 HRT 


0.074 


0.503 


-0.472 


-3.215 




-0.472 -3.215 


MC-EW 


0.057(1) 


0.539(3) 


-0.373(6) 


-3.55(2) 


-3 


.35(3) 






*(A«)„ 


- 1 (A 

NkT ^ 


A r " { ) v 










b (Aa), 


= n\t ( a ~ 


A»f), 
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Table 5. Critical parameters. Results based on the Carnahan-Starling hard sphere free energji^i for the reference system. 

0.0 OMF 1.13194 0.24913 



10~ 3 SCOZA 1.1319 0.2491 

0.1 SCOZA 1.1294 0.2495 

HRT 1.129 0.250 

OMF 1.13108 0.24915 

0.5 SCOZA 1.0718 0.2586 

HRT 1.073 0.260 

1.0 SCOZA 0.9264 0.2793 

HRT 0.925 0.279 

1.8 SCOZA 0.6527 0.3145 

HRT 0.650 0.314 

4.0 a SCOZA 0.173 0.3895 

HRT 0.175 0.394 

7.0 a SCOZA 0.0187 0.4575 

HRT 0.019 0.424 



a from ref. Iioj . Due to a different definition, the temperatures of ref. lipj must be multiplied by the factor ct* 2 /e a . 

^ r i a i , N sinh(ar) , . _. 

Q T = / dr47rr 2 T(r) — -. (A5) 

Jo ar 

As in the case of the Coulomb potential, the potential outside the spherical distribution 
of charge y(r) is still a Yukawa potential with the same screening parameter a but with 
a different charge Q T . Since sinh(x)/x > the effective charge Q T is larger than the bare 
charge r(0) of the distribution. Of course in the limit a —>■ one has Q T —>■ 1 and Gauss 
law is recovered. 

We now prove the theorem 

-^P + V(i2)>0 Vi2. (A6) 

For R > a it is obvious; for R < a we note that from the expression (1A5I) of Q T and from 
equation ( 1A2I) it follows that 

V T (R) - Q T y(R) = [ dr 4irr 2 T(r)^—Z(r) , (A7) 

J R OLVR 

where Z(r) = e~ ar smh(aR) — sinh(ar)e~ aR . Since Z(r) is a decreasing function of r and 
Z(R) = we have Z(r) < from which inequality fl A6[) follows (note that r(r) must be 
positive). For a sufficiently regular distribution r(r), the potential V T (R) is a bounded 
smooth function for R < a, in particular V T (R = 0) is finite in general (see examples at the 
end of the appendix). 
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Clearly the Fourier transform V T (k) = 47rr(fc)/(a 2 + k 2 ) has not a definite sign in general 
although K(0) = 47rf(0)ar 2 > 0. 

The interaction W T (1, 2) of two identical spherical distributions of Yukawa charges centred 
at points Ri and R2 is now given by (with obvious notations) 

W T (1, 2) = J dl' d2' r(l, 1') y(l', 2') r(2, 2') 

= J dl'K(M'Ml',2). (A8) 

Also in this case Gauss theorem generalises easily. Indeed, for the case where the two 
distributions do not overlap, i.e. R = |R 2 — Ri| > er, V T (1, 1') can be replaced by Q T y(l, V) 
in equation ( 1A8I) and thus 

W T (R) = Ql y(R) for R > a . (A9) 

Note that the Fourier transform W T (k) = r 2 ^) y(k) > is positive definite contrary to 
V T {k). 

In addition, we have the theorem 

- + y(R) > Vi? . (Aio) 

The proof is trivial : 



V(l,2) = J dl'r(l,l')V T (l',2) 

<Q T J dl' r(l, 1') y(l', 2) = Q T V T (1, 2) 



<Qiy(l,2). (All) 

The function W T (r) is in general bounded in the core and in particular W T (0) is finite. This 
charge smearing process provides an easy way to regularise the Yukawa potential y(r) in the 
core which preserves the positivity of the interaction. 

Finally, we give explicit expressions for Q T , V T for r < a, and the value W T (0) in simple 
cases relevant for Sect. IIVCI . Thus for a surface distribution 

Ta (r) = -^5(r-a) , (A12) 
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one finds 

sinh(aa) 



Q 



a 

* » exp(— aa) sinh(ar) „ 
V T (r) = — ^ — - for r < a 

a ar 

W T ,(0) = Qr, . (A13) 

For a volume distribution 

r p (r) = —Q(a-r), (A14) 



(6 step function) one has 



3 

Qt„ = _> o (ao r cosh(aa) — sinh(cKf)) 
(aa) 6 



j r / \ 3 / . _, smh(ar)\ _ 

K- fl (r) = 1 - (1 + aa) exp (-okj) — - for r < a 

p a 2 a 6 \ ar J 

W Tp (0) = -J^ (l-(l + aa)e W (-aa)Q Tp ) . (A15) 
a a 
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I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 

SCOZA 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

p* 

FIG. 1: Liquid-vapour coexistence curve of the HCY fluid for a* = 1.8 from SCOZA and HRT. 
The simulation results (open squares) are for the cubic system with EW sums (N = 1000). 
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0.9 - 



n — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — r 



SCOZA 

a- -a HRT 
□ MCC3 




FIG. 2: Liquid-vapour coexistence curve of the HCY fluid for a* = 1.0 from SCOZA and HRT. 
The simulation results (open squares) are for the cubic system with EW sums (N = 1000). 
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t — i — i — r 



SCOZA 




0.1 0.2 0.3 0.4 0.5 0.6 



FIG. 3: Liquid-vapour coexistence curve of the HCY fluid for a* = 0.5 from SCOZA, HRT and 
OMF. Open squares: simulation results for the cubic system with EW sums (N = 1000); open 
circles: hypersphere calculations (N = 2000). 
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i — i — i — r 




0.1 0.2 0.3 0.4 0.5 

p* 



FIG. 4: Liquid- vapour coexistence curve of the HCY fluid for a* =0.1 from SCOZA, HRT and 
OMF. Simulation results are for the cubic system with EW sums (open squares: iV = 1000, filled 
squares: N = 1872); and hypersphere calculations (open circles N = 2000). 
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~i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — r 
SCOZA 




-5 -4 -3 -2 -1 

U/NkT 



FIG. 5: Internal energy U/Nk^T along the coexistence curve for a* = 0.5. The symbols are as in 
Fig. 3. 
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